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We study the nuclear symmetry energy S{p) and related quantities of nuclear physics and nuclear 
astrophysics predicted generically by relativistic mean-field (RMF) and Skyrme-Hartree-Fock (SHF) 
models. We establish a simple prescription for preparing equivalent RMF and SHF parameteriza- 
tions starting from a minimal set of empirical constraints on symmetric nuclear matter, nuclear 
binding energy and charge radii, enforcing equivalence of their Lorenz effective masses, and then 
using the pure neutron matter (PNM) equation of state (EoS) obtained from ab-initio calculations 
to optimize the pure isovector parameters in the RMF and SHF models. We find the resulting 
RMF and SHF parameterizations give broadly consistent predictions of the symmetry energy J 
and its slope parameter L at saturation density within a tight range of < 2 MeV and < 6 MeV 
respectively, but that clear model dependence shows up in the predictions of higher-order symmetry 
energy parameters, leading to important differences in (a) the slope of the correlation between J 
and L from the confidence ellipse, (b) the isospin-dependent part of the incompressibility of nuclear 
matter Kt, (c) the symmetry energy at supra-saturation densities, and (d) the predicted neutron 
star radii. The model dependence can lead to about 1-2 km difference in predictions of the neutron 
star radius given identical predicted values of J, L and symmetric nuclear matter (SNM) saturation 
properties. Allowing the full freedom in the effective masses in both models leads to constraints of 
30 < J < 31.5 MeV, 35 < L < 60 MeV, -330 < < -216 MeV for the RMF model as a whole 
and 30 < J < 33 MeV, 28 < L < 65 MeV, -420 < < -325 MeV for the SHF model as a whole. 
Notably, given PNM constraints, these results place RMF and SHF models as a whole at odds with 
some constraints on Kt inferred from giant monopole resonance and neutron skin experimental 
results. 

PACS numbers: 21.65. Cd, 21.65.Mn, 26.60. Kp, 26.60.-c 



I. INTRODUCTION 



Highly isospin asymmetric nuclear matter is present in heavy nuclei far from the stability line and in the surface 
regions of nuclei exhibiting neutron skins and occurs during heavy-ion collisions and also in astrophysical systems 
such as neutron stars. Its energy relative to that of symmetric nuclear matter (SNM) is usefully characterized by the 
symmetry energy as a function of density S{p), and therefore the constraining of S{p), particularly at densities away 
from nuclear saturation density po, has been the subject of much recent experimental and theoretical activity [1-19]. 
Particularly, constraints extracted for the magnitude of the symmetry energy at saturation density J = S{pq) and its 
slope there L span the ranges J « 25 — 35 MeV and L « 25 — 115 MeV although the J constraints from mass models 
alone are much tighter, and most recent L constraints are placed in the range 30 — 80 MeV (see, e.g., and Fig. 
1 from [20j). Extraction of such constraints requires specifying a model for nucleon-nucleon interactions which tends 
to be equivalent to specifying a choice of functional form for S{p) in the RMF and SHF models. It is the influence of 
such choices on predictions of symmetry energy and simple, related, nuclear and astrophysical properties that will be 
the subject of this paper. 

Given the currently prohibitive complexity of computing the nucleon-nucleon (NN) interaction from the underlying 
theory of QCD and the challenges of solving the nuclear many-body problems, progress in nuclear many-body theory 
has developed along two lines: 

(1) The microscopic approach builds up the many body system from bare NN interactions together with a usually 
phenomenological description of the three-nucleon (3N) interaction; in-medium correlations are self-consistently in- 
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eluded via the many-body calculational technique. In the last few years, much theoretical progress in understanding 
neutron-rich systems has been accomplished through microscopic pure neutron matter (PNM) calculations. By study- 
ing the universal behavior of resonant Fermi gases with infinite scattering length, a significant constraint is achieved 
for the equation of state of dilute neutron matter |^ . These calculations have been extended to higher densities using 
the full power of Quantum Monte Carlo methods [521 112 ■ By studying the physics of chiral thrcc-nuclcon forces the 
EoS of PNM is obtained perturbatively up to nuclear saturation density [24] . Finally, the Auxiliary Field Diffusion 
Monte Carlo (AFDMC) technique, which takes into account the realistic nuclear Hamiltonian containing modern two- 
and three-body interactions of the Argonne potential and Urbana family of three-body nucleon forces, has been used 
to calculate the EoS of PNM up to and above saturation density [55H?f| . The outcome of the above investigations 
is a robust prediction for the EoS of PNM at low densities where only two body interactions are important, and a 
systematic investigation of the uncertainties in the EoS of PNM up to and beyond saturation density as a result of 
our present uncertainties in the three-neutron interaction, resulting in "theoretical error bars" in that density regime. 

(2) The second approach is to construct an effective interaction describing the in-medium nucleon-nucleon interac- 
tion, subject to most of the symmetries of the bare potential. The effective interaction is typically dependent on ^ 10 
parameters representing, for example, coupling constants, which are fit to experimental data sets from finite nuclei 
properties such as binding energies, charge radii, single particle energy spectra and spectra of collective excitations. 
One of the overriding goals of modern nuclear many-body theory is to derive an energy-density functional (EDF) |28j 
with clear physical connections to ah initio NN interactions and QCD. The widely used RMF [29l[30] and SHE [3TJ[32] 
models, with the latter thought of as a non-relativistic expansion of the former [33, 34^, are two typical phenomenolog- 
ical EDFs used in nuclear many-body theory. Both models have < 10 free parameters in their simplest forms. Recent 
surveys find about 240 parameter sets for the SHF model [35] and IDs of parameterizations of the simplest form of 
the RMF model, e.g. [3S1, although many are old parameter sets superseded by parameter sets fit to more recent, 
accurate data. We shall refer to the space inhabited by the free parameters as the model space, and the two classes of 
EDF (RMF and SHF) shall be referred to as the two models. Since the number of experimental observables is always 
larger than the number of free parameters, the problem of optimizing these EDFs is always overdetermined, and this 
results in a considerable degeneracy amongst parameter sets, and correlations between individual parameters when 
constrained by certain observables. Covariance analysis techniques [371 138| have been employed to study correlations 
between predicted observables from a particular EDF in its model space. Given a certain experimental constraint, 
this analysis method serves as a useful tool to optimize the parameters of RMF and SHF EDFs and systematically 
examine the correlations between various nuclear matter and neutron star properties jl6j. 

Most parameterizations of EDFs are obtained through fitting to predominantly nuclear experimental data sets, thus 
probing closely isospin symmetric systems; their predictions of symmetry energy behaviors thus vary widely, and many 
do not give PNM predictions consistent with microscopic calculations. Recent parameterizations which do take into 
account ab-initio PNM calculations tend to give behaviors of the symmetry energy at odds with some experimental 
constraints from giant monopole resonances and neutron skins [71 1391 HO] , constraints extracted, in part, using the 
same types of EDFs. This discrepancy also occurs in microscopic calculations [41]. This raises questions such as: can 
such discrepancies be resolved by choosing a different parameterization? Are the extracted constraints correct? Is 
there a fundamental problem with the particular EDFs used? A systematic study of parameterization-independent 
RMF and SHF model predictions of the behavior of J and L and related physical properties has yet to be undertaken. 

RMF and SHF models predict different functional forms for S{p) and thus one might expect generic differences 
in the values extracted from the same sets of data within each model, and conversely, generic differences in the 
predictions of properties of neutron-rich systems, even given the same values of J and L. This potentially makes 
combining the two types of functional in modeling experimental phenomena hazardous; it also means that extraction 
of, for example, constraints on J and L from nuclear experiment or astrophysical observation (e.g. the measurement 
of neutron star radii) comes with the caveat that such constraints are dependent on the model used for extraction. It 
is important to attempt to quantify what difference that model choice makes. 

The aim of this paper is to explore the generic predictions of properties of isospin-asymmetric nuclear matter from 
RMF and SHF EDFs simultaneously constrained by the best theoretical knowledge of the PNM EoS. Particularly, 
we will set model-generic best fit values and la confidence intervals on J, L, K-^ and neutron skin thicknesses arising 
from the optimization to the theoretical PNM EoS; we will examine whether the discrepancies between the predicted 
values of Kt and those extracted from experiment are endemic to the two EDFs as a whole, and we will explore the 
model dependence of the above results, the supra-saturation symmetry energy and neutron star properties arising 
from the different functional forms of the two models. We should note here that our aim is neither to establish new 
parameterizations of either model, nor to set absolute constraints on symmetry energy, but to explore as far as possible 
the generic constraints that can be placed by each model on neutron-rich systems once constrained by information 
from the PNM EoS. 

The manuscript is organized as follows. In Sec. [IT] we briefly review the two EDFs and the covariance analysis 
method used to optimize the two pure isovector parameters in the EDFs. Results are presented in Sec. [Hi] and in 
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Sec. |TV]we conclude. 

II. FORMALISM 

To ease discussions, we recall in this section the main formulas related to the symmetry energy, the RMF and SHF 
models, and the covariance analysis method used to examine the effect of the PNM constraints on the two models. 

A. Symmetry energy 

The binding energy per nucleon in neutron-rich nuclear matter can be written as 

E{p,a)^Ea{p) + S{p)a' + 0{a^), (1) 

where p is the baryon number density and a = (pn — P-p)/ P is the isospin asymmetry, with Pu{p-p) being the neutron 
(proton) number density. Around the saturation density po, the symmetry energy can be expressed as 

5(p) = J + Lx+^^symx' + 0(x') , (2) 

where x = (p — po) /3po: J is the value of the symmetry energy at saturation density, L is the slope parameter, and 
Ksym is the curvature parameter at saturation density given, respectively, by the following expressions: 

Ksym = 9po ( j _ • (4) 



P=PO 



The coefficients of the higher-order terms in Eq. ([T]) are generally much smaller than S{p), so it is usually a good 
approximation to write the energy per nucleon in PNM as -Epnm(p) ~ E!q{p) + S{p); however, in this work we calculate 
^'PNM using the full EoS. 



B. Relativistic mean-field model 



The commonly employed RMF model contains an isodoublet nucleon field {ip) interacting via the exchange of the 
scalar-isoscalar tr-meson (0), the vector-isoscalar w-meson (F^), the vector-isovector p-meson (b^), and the photon 
{A'^) [291 |30l |42] . The effective Lagrangian density for the model can be written as 

- Iv^'V^. + ImlV^V^ - ib^^ . b^, + ^-ml b'' • b^ - ^^^^''F^. - U{^, V^, b^) , (5) 
where V^j^, b^j^, and F^^ are the isoscalar, isovector, and electromagnetic field tensors, respectively: 

Vp, = d^V, - d,V,, , (6a) 
bp,. = d^h^ - d^hf, , (6b) 
Ff,, = d^A, - d,A^ . (6c) 

The nucleon mass M and meson masses m^, m^, and rrip may be treated as empirical parameters. The effective 
potential U{(j3,Vp_,h^) consists of non-linear meson interactions that simulates the complicated dynamics encoded in 
just few model parameters. In the present work we use the following form of the effective potential [43j : 

u{ct>,v^hn ^ ^{9s<f>r+^{9sci^r-^gt{v,vn' - A.^^^^ (7) 

This model is described by 7 interaction parameters: {gs, <?v, ffp, k, A, C, Av}. Note that power counting suggests that 
a consistent Lagrangian density should include all terms up to fourth order in the meson fields. However, the existing 
database of both laboratory and observational data appears to be accurately described by the the minimal set of 
parameters P3ti45j . Indeed, it was shown that ignoring a subset of model parameters that are of the same order in a 
power-counting scheme does not compromise the quality of the fit |421 146] . 
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C. Skyrme-Hartree-Fock model 

The standard form of the energy density obtained from the zero-range Skyrme interaction using the Hartree-Fock 
method can be written as \M\ 



2M 



-+to [(2 + x„) p2 _ (2x0 + 1) {pI + pI)] /4 



+ hp- [(2 + xg) p' - {2x3 + 1) {pI + pI)] /24 

+ [t2 {2X2 + 1) - tl (2X1 + 1)] {TnPn + TpPp) /8 + [ti (2 + Xi) + (2 + X2)] Tp/8 

+ [Ml (2 + xi) -t2{2 + X2)] {Vpf /32 - [3ti (2x1 + 1) + ^2 {'2x2 + 1)] 



(Vp„)V(Vpp)' /32 



Wo 



J- Vp+Jn- Vpn + Jp ■ Vpp /2 + {h - t2) [Jl + Jp] /16 - (tiXi + t2X2) / 1& 



(8) 



Here p^, r^, and Jq {q = p, n) are, respectively, the number, kinetic, and spin-current densities, and p, r, 
and J are the corresponding total densities. The SHF model is expressed in terms of 9 Skyrme parameters: 
{to,ti,t2,t3, xo,xi,X2, X3, a} and the spin-orbit coupling constant Wq which is taken as 133.3 MeVfm^ [10^ in the 
present work. 

D. Covariance analysis method 

Here we very briefly discuss the covariance analysis method used in the present work. For more details, we refer the 
readers to Refs. [371 [351 SB]- Given a set of N experimental observables On^^^ that are determined with an accuracy 
of AOn, one can minimize the quality measure x^- 



X'(P) 



N 

E 

n=l 



0i^'^^(p)-drp)' 



(9) 



Here each of the N observables is computed within the given model ^^^(p) as a function of the F model parameters 
p = (pi, . . . ,pf)- a set of optimal parameters po are determined via a least square fit to the quality measure. For 
our set of 'experimental' observables On^^'' in the input we choose the theoretical calculations of the energy per 
neutron Epjssm in the density range of 0.04 < p < 0.16 fm"^ [24l [25l |49] . Although the AFDMC calculations have 
been extended up to several times the saturation density |27| . the extension of the calculations of the chiral three- 
nucleon forces to higher densities using piecewise polytropes [5D] shows that the uncertainties in the EoS could be 
very large when all of these models are employed. Therefore we rely on the PNM calculations that are obtained up to 
saturation density only. Moreover, the symmetry energy coefficients are only sensitive to the equation of state around 
the saturation density. 

Once the optimal parameter set po is found through the x^-minimization, one can then compute and diagonalize 
the symmetric matrix of second derivatives. All the information about the behavior of the function around the 
minimum is contained in this matrix. That is. 



(10) 



where 



(P - Po) 
(Po)» 



(11) 



are scaled dimensionless variables, TW = ADA^ ^ and £^~A^ni are dimensionless variables in a rotated basis. Here A 
is the orthogonal matrix whose columns are composed of the normalized eigenvectors and ■D = diag(Ai, . . . , Xp) is the 
diagonal matrix of eigenvalues. The meaningful theoretical uncertainties can be obtained by computing the statistical 
covariance of two observables A and B which can be written as follows: 



cov(A, B) = ^ 



dA - dB 

{M 



dx. 



dxj 



(12) 
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The variance cr^(A) of a given observable A is then simply given by a'^{A)=cov{A, A). Finally, the covariance ellipses 
between two observables A and B can be plotted by diagonalizing the 2x2 covariance matrix: 

^ _ f cov{A,A) cov{A,B)\ , . 

cov(B, A) cov(B, B) J ^^•^> 

Then the eigenvalues of this matrix represent the semi-major and semi-minor axes of the covariance ellipse, while the 
eigenvectors provide the orientation of the ellipse. 



III. RESULTS 



Following the idea of building relations between values of model parameters and macroscopic nuclear quantities [lOj , 
one finds that by changing the two solely isovector parameters, which are gp and Av in the RMF model [5T], and xq 
and ^3 in the SHF model, only the symmetry energy S{p) is modified while properties of SNM such as saturation 
density pq, binding energy per nucleon at saturation density Eq, incompressibility coefficient at saturation density 
Kq^ and effective mass M* at saturation all remain unchanged. Thus, in the following we optimize the two isovector 
parameters [i^ = 2 in Eqs. (10 1 and (12|] with respect to the available range of PNM EoSs to constrain the values of 
J and L at saturation density by employing the covariance analysis method discussed above. 



A. Reference models 



As representative RMF parameterizations, we choose the accurately-calibrated NL3* [52] and the recent lU-FSU [55] 
parametrizations. The lU-FSU is the recent parameterization that was validated against experimental, observational, 
and theoretical data, while the accurately-calibrated NL3* parameterization gives a much stiffer EoS of SNM (larger 
value of Kq and smaller value of C parameter) and a stiff symmetry energy (larger values of symmetry energy J and 
slope L) and therefore offers a suitable contrast to lU-FSU. 

To compare the RMF and SHF models on the same footing, we create two Skyrme parameterizations which give 
the same properties of nuclear matter at saturation as the two RMF parametrizations, herein referred to as SkNL3* 
and SklU-FSU forces, through the method of writing the Skyrme parameters as functions of macroscopic nuclear 
quantities [101 154] . Note that these two new Skyrme parameterizations are intended only to serve as references in this 
study. 

Several definitions of the nucleon effective mass exist in the literature [55]. In the RMF model the Dirac effective 
mass is defined through the scalar part of the nucleon self-energy in the Dirac equation: 

Mu,, = M, + , (14) 

where the nucleon self-energy is given as Sf^ = Ep = —gs<t> in the RMF model considered in this work. It has 
been well documented that there is a strong correlation between the Dirac effective nucleon mass at saturation 
density and the strength of the spin-orbit force in nuclei [SO] [33l [56l Indeed, one of the most compelling 
features of RMF models is the reproduction of the spin-orbit splittings in finite nuclei. This occurs when the velocity 
dependence of the equivalent central potential that leads to saturation arises primarily due to a reduced nucleon 
effective mass It is shown that models with effective masses outside the range 0.58 < M^/M < 0.64 will not be 
able to reproduce empirical spin-orbit couplings 1581 . when no tensor couplings are taken into account. On the other 
hand, the non-relativistic effective mass parameterizes the momentum dependence of the single particle potential, 
which is the result of a quadratic parameterization of the single particle spectrum. A recent study |35j puts a bound 
of 0.69 < M* /M < 1.0 for the non-relativistic effective masses. It has been argued [5^] that the so-called Lorentz 
mass Ml should be compared with the non-relativistic effective mass extracted from analyses carried out in the 
framework of nonrelativistic optical and shell models. For consistency, we choose the effective mass in the Skyrme 
parameterizations to be equal to the Lorenz mass in the RMF parameterizations. Since the RMF model we use in 
this work gives the same isoscalar and isovector effective masses, we set them equal in the reference SHE model too. 

Finally, the isoscalar parameters of the two reference Skyrme forces are then re-adjusted to fit the binding energy 
and charge radius of ^osp^-, -^y adjusting only the saturation density po and the binding energy Eq of SNM. As shown 
in Fig. [TJ these models predict the charge radii and binding energies of other doubly closed-shell nuclides within 
1-2% accuracy. We note that these finite nuclei properties are obtained by solving the Dirac equation for the RMF 
model and the Schrodinger equation for the SHF model. The bulk nuclear matter observables predicted by these 
reference models are given in Table |lj In terms of the predicted values of isoscalar and isovector bulk observables, 
both corresponding RMF and SHF models are therefore almost equivalent. 
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FIG. 1: (Color online) Relative deviation of the binding energies (a) and charge radii (b) of closed-shell nuclei from the reference 
models discussed in the text compared with the experimental data (with superscript 'exp') from Refs. |60II61| . Filled symbols 
are from the original parameterizations and empty symbols are from the PNM modified parameterizations. 





Po (fm-^) 


Eo (MeV) 


Ko (MeV) 


Mi (M) 


Mt (M) 


Ms* (M) 


(M) 


J (MeV) 


L (MeV) 


iCeym (McV) 


i?skin (fm) 


NL3* 


0.1500 


-16.32 


258.49 


0.594 


0.671 






38.7 


122.7 


105.7 


0.29 


SkNL3* 


0.1527 


-15.76 


258.49 






0.671 


0.671 


38.7 


122.7 


62.7 


0.27 


lU-FSU 


0.1546 


-16.40 


231.33 


0.609 


0.687 






31.3 


47.2 


28.5 


0.16 


SklU-FSU 


0.1575 


-15.70 


231.33 






0.687 


0.687 


31.3 


47.2 


-132.0 


0.16 



TABLE I: Macroscopic quantities from four reference parameterizations. They are binding energy per nucleon Eq and incom- 
pressibility Kq of SNM, Dirac (Mp) and Lorentz (M£) effective mass from the RMF model, non-relativistic isoscalar (Ms) and 
isovector (My) effective mass from the SHF model, the symmetry energy J, its slope parameter L and curvature parameter 
Ksyui at saturation density, and the resulting neutron skin thickness RaMn of ^''^Pb. 



The energy per neutron ii'pNM predictions at sub-saturation densities for our reference models are plotted on the 
left panel (a) of Fig. [2] compared to the results obtained by various microscopic approaches. One can see that even 
among our four parameterizations there is wide variance in the EoS of PNM at all densities, and little agreement 
with those microscopic calculations. The very wide range of predictions of the symmetry energy parameters and the 
corresponding widespread predictions for the neutron skins of nuclei inherent in these parameterizations are seen in 
Table H 



B. Symmetry Energy Coefficients 

Having established our reference models, we next minimize the with respect to the PNM constraints [211 HUSH] 
in the density range of 0.04 < p < 0.16 fm~^ by adjusting two isolated (solely isovector) parameters. This leads 
to optimized values of the model parameters and thus the density dependence of symmetry energy up to saturation 
density once the EoS of SNM is fixed. All isoscalar parameters remain unchanged and there is very little change in 
the prediction of binding energies and charge radii as is shown in Fig. [l] 

As can be seen in panel (b) of Fig. [2j we obtain the EoS of PNM for a given RMF or SHF parameterization that 
best fits within the band of microscopic PNM calculations at the minimum of the x^-function. The resulting RMF and 
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FIG. 2: (Color online) Comparing the PNM EoS from four reference parameterizations with the AFDMC EoS in the AV8'+UIX 
Hamiltonian [27], the variational APR EoS [49], the low-density band from the constraints of resonant Fermi gases [21], and 
the high-density band from the chiral effective field theory calculations with 3- neutron forces [55 , before (a) and after (b) PNM 
optimization. 





^0.1 


So.i 




J 




L 


'-sym 


-^sym 






-'•'skin 


-Rs]<;in 


NL3* 


25.7 


24.9 ± 0.4 


38.7 


30.7 ± 0.7 


122.7 


50.3 ± 1.8 


105.7 


39.2 ± 17.8 


-684.4 


-284.6 ± 29.4 


0.29 


0.18 ± 0.01 


SkNLS* 


25.0 


24.5 ± 0.3 


38.7 


31.0 ± 0.9 


122.7 


46.4 ± 6.4 


62.7 


-156.0 ± 16.6 


-529.4 


-380.0 ± 15.2 


0.27 


0.16 ± 0.01 


lU-FSU 


25.7 


24.9 ± 0.4 


31.3 


31.4 ± 0.7 


47.2 


52.9 ± 2.0 


28.5 


-6.8 ± 12.9 


-195.3 


-257.6 ± 22.3 


0.16 


0.18 ± 0.01 


SklU-FSU 


24.4 


24.4 ± 0.3 


31.3 


31.4 ± 0.9 


47.2 


48.0 ± 6.2 


-132.0 


-130.2 ± 13.3 


-343.9 


-345.6 ± 15.3 


0.16 


0.16 ± 0.01 



TABLE II: Isovector observables and associated la error bars from four reference parameterizations before (with superscript '0') 
and after (without superscript '0') the PNM constraints are applied. Values are shown for the symmetry energy at p = 0.1 fm~^ 
So.i and at saturation density J, slope parameter L, curvature parameter -ft'sym, isospin-dependent part of incompressibility 
Kt, and the neutron skin thickness i?akin of '^'^^Pb. All the quantities are in MeV apart from -Rskin which is in fm. 



SHF models predict very similar symmetry energies J, while the RMF model predicts a consistently higher central 
value for L by about 4-5 MeV than the SHF model as shown in Table [n] 

The la errors on these two isolated parameters can be translated into equivalent errors on J and L. The errors in 
J are less than ±1 MeV for all the parameterizations. The RMF model gives a relatively small error in L of around 
±2 MeV, while the SHF model gives a much larger error around ±6 MeV. Table |Tl] appears to indicate that within 
the la errors, both models are consistent in their predicted values of J and L. However, in Fig. |3]we plot a la joint 
confidence regions in the J-L plane for both RMF and SHF models, thus showing that in fact the two models predict 
non-overlapping regions in J-L space. Both models show a positive correlation between J and L, but with differing 
slopes. For example, for lU-FSU and SklU-FSU parametrizations the relations are approximately 



L = (2.4 J- 23) MeV, (RMF) 
L = (6.0 J -140) MeV, (SHF) 



(15) 
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within the constraints of J and L shown in Table ITll 

The origin of this difference hes mainly in the values of the higher-order symmetry energy parameters that are 
predicted upon optimization. There is a strong model dependency in the prediction for the curvature parameter of the 
symmetry energy fCsym (see Table [Tl|) . For example, after the PNM optimization lU-FSU predicts i^sym = — 6.8± 12.9 
MeV, while its Skyrme-like version predicts a smaller value of Kgym = —130.2 ± 13.3 MeV. When wc plot the Ict 
joint confidence regions in the Kgym-L plane for both RMF and SHF models (see the left panel (a) of Fig. |4]) further 
differences can be seen: there is, generically, a negative correlation between the slope of the symmetry energy and 
Ksyn\ in the RMF model, while this correlation is positive in the case of the SHF model. Only at a sub-saturation 
density of p = 0.1 fm^'^ do the two models have similar values of ii'sym (p = 0.1 fm^"^) (see the right panel (b) of Fig. 
|4]), although the correlations are still opposite. We emphasize that these qualitative features emerge whatever the 
starting parameterization of the RMF or SHF model used. 




29 30 31 32 33 29 30 31 32 33 

J (MeV) 



FIG. 3: (Color online) la joint confidence regions for the symmetry energy J and its slope parameter L at saturation density 
for the RMF and SHF models. 



It is widely accepted that the Giant Monopole Resonance (GMR) provides the cleanest and most direct route to 
the nuclear incompressibility around normal density j62) . It has been also proposed that GMR energies of finite nuclei 
as well as the nuclear matter incompressibility should be computed within the same theoretical framework [631 164) . 
The expression for the incompressibility of neutron-rich matter at saturation density is given by [SI] : 

i^sat(a) ^Ko + Kro" + 0{a^) , (16) 

where the coefficient of is 

Kr^Ksy^-6L^^L (17) 

with Qo being the skewness of SNM [54j. Although both RMF and SHF models used in this work share the same 
value of Kq, their predictions of -K^sat ai'e different due to the difference in K^^ which in turn is mainly due to the 



difference in -fsTsym- In Table HI we provide the values of Ksat for different values of isospin asymmetry. Due to the 
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FIG. 4: (Color online) la joint confidence regions for the slope parameter L and curvature parameter A"sym of the symmetry 
energy at saturation density (a) and at p — O.l fm~^ (b) from the lU-FSU and SklU-FSU parameterizations. 



small values of isospin asymmetry in finite nuclei, the difference of the incompressibility for different models is in fact 
small. Comparing with the constraint of —760 < Kr < —372 MeV in Ref. [35] extracted directly from the GMR data, 
both RMF and SHF models predict marginally consistent or slightly higher values of Kt after the PNM optimization 
as shown in Table |ll) suggesting that both RMF and SHF models have difficulty in simultaneously predicting GMR 
properties consistent with experiment and the PNM EoS consistent with our best theoretical calculations. 





K,^t{a^O) (MeV) 


A^sat (a 


= 0.111) (MeV) 


-ft'sat {a 


= 0.212) (MeV) 


NL3* 


258.5 




255.0 




245.8 


SkNL3* 


258.5 




253.8 




241.5 


lU-FSU 


231.3 




228.2 




219.8 


SklU-FSU 


231.3 




227.1 




215.9 



TABLE III: Incompressibility of neutron-rich matter with different isospin asymmetries a = (SNM), 0.111 (^"Zr) 
(^"^Pb) at saturation density from the four parameterizations after the PNM optimization. 



and 0.212 



Different values of the bulk properties of SNM will affect the PNM constraints on the symmetry energy. For example, 
the saturation density po, the binding energy at saturation Eq, and the incompressibility coefficient at saturation Kq 
will affect the EoS of SNM and thus modify slightly the optimized symmetry energy from a fixed set of PNM EoS 
constraints. The effective mass M* dominates these uncertainties in the results of the PNM optimization. As can be 
seen from the expression for the symmetry energy in the RMF model ^65j : 

M* affects the kinetic contribution to the symmetry energy while adjusting gp and Av only modifies the potential 
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J (MeV) 


L (MeV) 


(MeV) 


RMF 
SHF 


30.2 — 31.4 
30.1 — 33.2 


36.1 — 59.3 
28.5 — 64.4 


-329.7 — -215.7 
-418.8 — -235.3 



TABLE IV: Predicted ranges for symmetry energy parameters within RMF and SHF models with their pure isovector parameters 
optimized to PNM and taking into account all remaining variation from parameterizations constructed since 1995. 

contribution to the symmetry energy. We find that increasing the effective mass at saturation by ~ 10% decreases 
the optimized value of the slope of the symmetry energy at saturation density L by ~ 10 MeV. The isovector effective 
mass, here set equal to the isoscalar effective mass in the SHF model to be consistent with the RMF models, affects 
the value of L obtained in the PNM optimization by the same order of magnitude, but in the opposite direction. The 
curvature of the symmetry energy if gym is changed by a much smaller relative amount. Therefore the l-cr confidence 
ellipses change their positions in the J-L plane as the SNM properties are varied, but they retain very similar values 
of their slopes, and the RMF and SHF confidence ellipses maintain their relative positions. Similarly, the Ksym-L 
confidence ellipses change their L-position upon variation of SNM properties, but retain their i^sym values and relative 
orientation and spacing. 

In order to get a better idea of the overall range of predictions for J, L and AV taking the additional model 
parameters into account, we take 11 RMF parameterizations and 73 SHF parameterizations from the literature that 
have been created since 1995 |35l 136 ] . We optimize the pure isovector parameters of each parameterization to the 
PNM results and examine the resultant constraints; these are displayed in Table |IV] 

C. Implications for predictions of neutron skin thicknesses and neutron star radii 

Measurements of the neutron skin thicknesses of various nuclides using strong interaction probes [H and, 
recently, weak interaction probes |75l 176] in the PREX experiment, are an important tool to probe the density 
dependence of the symmetry energy due to the very close correlation of L to the size of the neutron skin in neutron- 
rich nuclides [71 HTl 177^179) . Since our optimized RMF and SHF models give nearly matching ranges of L, we expect 
the neutron skin predictions to be similar. 

In Fig. [5] we compare predictions of neutron skin thicknesses from the lU-FSU and SklU-FSU parameterizations 
to the currently existing data on the neutron skin thickness of Tin isotopes [H [671 17lll74j . As expected, both the 
post-optimization lU-FSU and SklU-FSU models agree well with the experimental data, with the RMF model giving 
a systematically slightly higher value than the SHF model in all but the lightest isotopes calculated. Thus consistency 
with our best knowledge of the PNM EoS can be achieved simultaneously with consistency of neutron skin predictions 
with current experimental data within the RMF and SHF models. 

The lU-FSU parameterization predicts i?skin = 0.18 ± 0.01 fm for ^''^Pb, while SklU-FSU predicts a slightly lower 
value of i?skin = 0.16 ± 0.01 fm (Table |ll]). The smaller value of i?skin for SklU-FSU is primarily due to model 
dependence, which leads to a smaller value of optimized L from the PNM constraints. The current PREX obtained 
value for the neutron skin thickness of lead is i?skin = 0.33^° ;[g fm [75]. If the new PREX experiment reduces the 
error bars without moving the central value for the neutron skin, almost all current models of the nuclear structure 
would need to be modified. Also, this would appear to call for a significant modification of the PNM microscopic 
calculations. 

Finally, we examine how the different symmetry energy characteristics of RMF and SHF models are manifest in 
neutron star radius predictions. Using our four post-optimization parameterizations, we apply the EoS of /3-stable 
and charge neutral neutron star matter composed of neutrons, protons, electrons, and muons throughout the core of 
the star. For the very low density outer crust we use the BPS equation of state [80]. The equation of state of the 
inner crust is approximated by the polytropic equation of state of the form P = AS'^'^ + B [5T] , where A and B are 
determined to match the EoS at the boundaries of the inner crust. Using our equations of state, we integrate the 
general relativistic equation for hydrostatic equilibrium (the Tolman-Oppenheimer-Volkoff equation) from the center 
to the surface of the star. 

The reference RMF and SHF parameterizations before the PNM optimization predict a wide range of results for 
low mass neutron star radii as shown in the left panel (a) of Fig. |6j In particular, for a 1.0 solar mass neutron star the 
difference in the predictions of radii for the NL3* and lU-FSU equations of state is equal to Ai?i.o ~ 2.8 km. There 
is a similar difference between the original SkNL3* and SklU-FSU equation of state predictions, i.e., Ai?i.o ~ 2.5 
km. This can be mainly attributed to the density dependence of the symmetry energy, which is quite different in the 
two parameterizations. Once calibrated to the PNM results, this difference almost vanishes within the same model 
as shown in the right panel (b) of Fig. [61 i.e., both RMF and SHF parameterizations now match each other more 
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FIG. 5: (Color online) Comparing the predictions of the neutron skin thickness for Sn isotopes from the lU-FSU and SklU-FSU 
models after the PNM optimization with those from different experimental methods. 



closely (excepting the differences at high masses between the RMF models, a result of a stifFer EoS of SNM in NL3* 
parameterizations at several times saturation density due to the C parameter). Although both NL3* and lU-FSU 
parameterizations in a given RMF or SHF model predict similar radii, there is a clear difference between the RMF 
and the SHF predictions as a whole. In the case of lU-FSU and SklU-FSU we have almost a 1 km difference for 
the radius of a canonical neutron star. This discrepancy is even larger in the case of NL3*, which is about ^1.8 km. 
Thus, there is a strong model dependence when the two models are applied to neutron star structure calculations 
after the same PNM optimization. 

The above model dependence actually comes from different density dependence of symmetry energy at supra- 
saturation densities, flagged by the model dependent difference in predictions of the curvature of the symmetry 
energy at saturation density i^sym- In Fig. IT] we plot the density dependence of the symmetry energy for the RMF 
and SHF models under consideration after the PNM optimization. Note that the symmetry energy is almost the 
same in all the models up to ^ 1.5/9o saturation density. However, the symmetry energy in the RMF functional is a 
monotonic increasing function of density, while the SHF functional tends to give a decreasing symmetry energy with 
increasing density at higher densities. Again, this property is generic once the model has been optimized to PNM 
EoS. The reason for this difference is manifest in the functional forms of the symmetry energy given as: 

5RMF(p) = A(p)p2/3 + B(p)p, (19) 

SsMp) = ap"^ -bp- cp^'^ - dp-+' , (20) 



2 / 2 \ ^/"^ 

where A{p) and B{p) are positive-valued functions of density [see Eq. (18)], a = ^ ( j and b, c, d are constants 
that depend on Skyrme parameters only. The symmetry energy in the RMF model is always positive as given in Eq. 



, while certain terms of the symmetry energy in the SHF model can become negative at higher densities [see Eq. 



Recently, it was shown that currently available neutron star mass and radius measurements provide significant 
constraints on the EoS of PNM all the way up to several saturation densities [82] . While this is true, we also show 
that the low-density PNM constraints alone result in a pronounced model dependency of radius predictions, as different 
masses and radii can be obtained with the similar saturation properties constrained by the low-density PNM EoS. 
Although our PNM optimization tightly constrains the symmetry energy up to a little above the saturation density. 
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FIG. 6: (Color online) Mass-Radius relation of neutron stars calculated from the four parameterizations before (a) and after 
(b) the PNM optimization. 



in order to understand its behavior at higher densities, which is also important in determining neutron star radii, one 
must rely on the heavy- ion collision experiments |83H85j and neutron star observations [861 187j . 



IV. CONCLUSIONS 



Using parameterizations of RMF and SHF energy-density functionals prepared to give equally good fits to ground 
state properties of doubly magic nuclei and identical symmetric nuclear matter properties, and are fit to state-of-the- 
art ab initio theoretical calculations of PNM up to saturation density, we have conducted a systematic examination 
of the resultant predictions from both models of the symmetry energy as a function of density and some important 
terrestrial nuclear and neutron star observables sensitive to S{p). 

We show that such RMF and SHF models result in very similar predictions for the symmetry energy J and its 
slope parameter L at saturation density from both models so long as the isoscalar effective mass from the SHF model 
is chosen to be equal to the Lorenz effective mass from the RMF model, which is tightly constrained around « Q.7AI. 
Both models then give J « 31.0 ± 1 MeV. The SHF parameterizations give values around 46 — 49 ±6 MeV and 
the RMF parameterizations 50 — 53 ±2 MeV for L. Resulting predictions of neutron skin thicknesses -Rskin for Sn 
isotopes and ^"^Pb therefore agree closely and are consistent with the available experimental data. 

When the la error bounds are plotted as ellipses in the J-L plane, a positively-correlated relationship between J 
and L is observed for both models. However, different slopes are obtained from the RMF and SHF models, and the 
two ellipses have no overlapping area in the plane. This model dependence comes from the different values of i^sym 
and higher-order symmetry energy parameters; i.e. from the different functional form of the symmetry energy implicit 
in the models. Although the PNM constraints lead to broadly similar behaviors of the symmetry energy as a function 
of density up to w 1.5po, they deviate significantly at higher densities due to the differences in the functional form of 
the symmetry energy. With the same PNM constraints up to the saturation density, the RHF model tends to predict 
a rising symmetry energy at higher densities, whereas the SHF model predicts a symmetry energy that may decrease 
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FIG. 7: (Color online) Density dependence of symmetry energy from the four parameterizations after the PNM optimization. 

with density at higher densities, and thus leading to the uncertainty of up to ~ 2 km in neutron star radii. Care must 
therefore be taken in extracting constraints on the symmetry energy, particularly on J and L, from inferred neutron 
star radii within one particular model. 

The absolute values of the predictions are found to be mainly sensitive to the effective mass, with increases (de- 
creases) of ~ 0.1 M leading to decreases (increases) of L by ~ 10 MeV. We confirm this systematic analysis by 
analyzing the predictions from 11 RMF and 73 SHF parameterizations constructed since 1995, finding overall ranges 
taking into account remaining freedom in the parameter values, of of 30 < J < 31.5 MeV, 35 < L < 60 MeV, 
-330 < Kr < -216 MeV for RMF models and 30 < J < 33 MeV, 28 < L < 65 MeV, -420 < Kr < -325 for SHF 
models. 

Notably, some recent constraints inferred from experimental data on giant monopole resonances of Sn and Cd 
isotopes [40l l88l [89] and on neutron skins 7 place Kr in the overall range —650 < Kr < —375 MeV. It has been 
pointed out that these results are inconsistent with many individual Skyrme parameterizations and microscopic nuclear 
matter calculations [33 [33] ; our results generalize these points to demonstrate that these particular Kr constraints 
are inconsistent with the RMF model as a whole and only marginally consistent with SHF models as a whole, within 
1(7 confidence intervals resulting from optimization to PNM calculations. Thus, either the density dependence of 
RMF and SHF models is insufficient to simultaneously describe PNM within current bounds and GMR/neutron skin 
experimental data, or there are overlooked problems with the extraction of the Kr constraints in the above works and 
the error bounds are underestimated, as has been suggested [31] . Note that the Kr ranges we extract from both models 
are consistent with another Kr constraint extracted from isospin diffusion in heavy ion collisions —490 < Kr < —250 
MeV [54; • More work needs to be done to check these hypotheses while taking the dependency of functional forms 
(e.g. [nO]) into consideration. 
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